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ABSTRACT: This article reviews a new and general theory of nonuniform fluids that naturally incorporates 
molecular scale information into the classical van der Waals theory of slowly varying interfaces. The method optimally 
combines two standard approximations, molecular (mean) field theory to describe interface formation and linear 
response (or Gaussian fluctuation) theory to describe local structure. Accurate results have been found in many 
different applications in nonuniform simple fluids and these ideas may have important implications for the theory of 
hydrophobic interactions in water. 
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This article reviews recent progress we and our coworkers have made in developing a new and 
general theory of nonuniform fluids (1-7), based on a reexamination of the ideas that lead to the 
classic van der Waals (VDW) theory (8, 9) of the liquid-vapor interface. The VDW interface 
theory, developed twenty years after the VDW equation of state for the uniform fluid, is equally 
far-reaching and has a great many virtues that merit further consideration in the light of modern 
developments in statistical mechanics. It is physically motivated, treating separately the excluded 
volume effects associated with the short-ranged and harshly repulsive intermolecular forces and the 
averaged effects of the longer ranged attractive interactions. Both thermodynamic and structural 
features are connected together naturally in an elegant and self-consistent approach. 

In principle the VDW theory can be applied to a very general problem: the structure and 
thermodynamics of a fluid in the presence of a general external field. However, the usual theory 
makes a crucial assumption that the fluid density in the presence of the field is in some sense slowly 
varying. While this assumption seems appropriate for the liquid-vapor interface in zero field where 
the VDW theory had its original and spectacularly successful application (9), it fails badly when 
applied to the more general and often rapidly varying fields associated with a number of problems 
of current interest. By trying to understand precisely where and why the classical theory failed, we 
have been able to develop a new perspective that allows us to address these more general problems. 

For example, the field can describe the interaction of fixed solutes with a solvent liquid. The 
solvation free energy is related to the changes in free energy as the interaction field is "turned 
on" from zero to full strength (10). A significant part of the solvation free energy arises from 
the required expulsion of solvent molecules from the region occupied by the solute. This involves 
very strong "excluded volume" interactions that can significantly perturb the density around the 
solutes. These considerations also play a major role in the theory of hydrophobic interactions (4, 
11-29), and we will later review and clarify work (4) by Lum, Chandler, and Weeks (LCW) based 
on this perspective. (As emphasized by Lawrence Pratt in his review of hydrophobic interactions 
in this volume (30), the general theory of hydrophobic interactions is incredibly complicated, and 
we will only deal with certain important but limited aspects here.) Similar issues arise in trying to 
understand effective interactions between solutes arising from depletion forces (31). 

An external field representing one or more fixed solvent particles leads to theories for multi- 
particle correlations functions describing the molecular scale structure of the bulk solvent liquid 
(32). Nonuniform fluids confined by walls, slits, pores, etc., can be described using appropriate 
external fields. Particular fields can enhance and alter features seen in the local structure of the 
uniform fluid. Confining fields can also produce shifts and rounding of bulk liquid- vapor and critical 
point phase transitions and can change the effective dimensionality of bulk fluid correlations; e.g., 
from three to two dimensional behavior in sufficiently narrow slits. In addition the field can induce 
new density correlations associated with interface formation, leading to wetting or drying transitions 
and related phenomena like capillary condensation (33-42). Thus, in an example we will consider 
in some detail later, a vapor-like drying region of lower density can form near a hard wall in a 
Lennard-Jones (LJ) fluid over a range of thermodynamic states near coexistence. We refer to the 
recent review by Gelb et al (33) for a detailed discussion of many of these possibilities. 
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This possible mixing of local structural components along with interfacial components in the 
density response to an external field leads to the great variety of behavior seen in nonuniform 
fluids. The classical VDW theory gives a qualitatively accurate description of the smooth interfacial 
components but fails for the local structure features that also must be taken into account for even 
a qualitatively accurate description of problems like those discussed above. 

Of course, it has long been recognized that the VDW assumption of a slowly varying interface is 
at fault. But that realization alone does not give us much insight into how to improve things. Most 
modern attempts (43-60) to address these problems have used density functional theory (DFT). 
Here one tries to express the free energy as a (generally nonlocal) functional of the density, and 
the VDW theory emerges when a local density approximation, appropriate for a slowly varying 
interface, is made. But what is that functional when the local density is rapidly varying and even 
gradient corrections are inadequate? Most workers have used a weighted density approximation in 
which one averages in some way over the rapidly varying local structure, and this approach has 
had great success in some applications. However, there is no "theory of theories" (61) for how to 
choose suitable weighting functions, and a host of different and often highly formal schemes have 
been proposed. These complications stand in strong contrast to the simplicity and physical appeal 
of the original VDW theory. We will discuss DFT further in Section [| 

We show here that there is another and very simple way to derive the VDW interface equation 
directly without first approximating the free energy. This allows one to think about the "local 
density" or "slowly varying profile" approximation in a new way where the theory itself suggests 
what is needed to correct it. Our work can be viewed as simply implementing these corrections to 
the classical VDW interface equation and it seems appropriate to refer to it for the purposes of this 
review as a molecular scale van der Waals (MVDW) theory. 

More precisely, we show that there are two key approximations in this interpretation of the VDW 
theory: i) the introduction of an effective single particle potential or "molecular field" to describe 
the locally averaged effects of the attractive interactions in the nonuniform fluid and ii) the use of 
a hydrostatic approximation to determine the density response to the effective field that takes into 
account only the local value of the field. The latter approximation is accurate only for very slowly 
varying fields, where it reduces to the local density approximation, and this is what limits the 
utility of the VDW theory in most of the applications discussed above. But from this perspective, 
it is easy to see how the theory can be corrected by considering nonlocal effects from the effective 
field. Our new MVDW theory does this using what is probably the simplest possible method, 
linear response theory (10, 62), which we argue (5, 6) is especially accurate when used to correct 
the hydrostatic approximation. 

In the following, we will present this interpretation of the VDW theory, first for the nonuniform 
LJ fluid. Then we will describe the corrections and generalizations that lead to the MVDW theory, 
and review results for several applications to nonuniform LJ and hard sphere fluids, and to water. 
Not all the simplifications that arise from this viewpoint were realized originally, so the present 
review can serve as a simpler and more concise guide to the theory. In particular, we have tried to 
clarify certain aspects of the LCW theory (4) for hydrophobic interactions. A different view of the 
LCW theory that may go beyond the MVDW picture considered originally and discussed here is 
given in (63). 
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We first consider the LJ fluid, where the physics is particularly clear. One reason why the LJ 
fluid is aptly called a simple liquid is that its local structure at typical liquid densities can be well 
approximated by an even simpler model, the hard sphere fluid. As pointed out by Widom (64, 
65), in most typical configurations of the uniform LJ fluid the vector sum of the longer-ranged 
attractive forces on a given molecule from pairs of oppositely situated neighbors tends to cancel. 
This leaves only the excluded volume correlations induced by the harshly repulsive molecular cores, 
well described by hard spheres of an appropriate size (10, 66). The hard sphere model is thus of 
fundamental importance in the theory of liquids as the simplest model that can give a realistic 
description of structural correlations arising from excluded volume effects. 

As emphasized by VDW himself (9), this cancellation argument must fail for nonuniform fluids, 
and there will exist "unbalanced" attractive forces, whose effect on the structure can also be quite 
significant, leading in some cases to interface formation. In general there is a competition between 
these two sources of structural correlations that a proper theory of nonuniform fluids should account 
for (1-7). This physical picture will help us understand what is missing in the classical VDW theory 
and how it might be improved. The insights gained from a careful study of this simple system may 
suggest physically well-grounded approximations that could be applied more generally. We believe 
this is the case for the MVDW theory. 

We start by defining the interactions in the nonuniform LJ fluid. Particles interact with a known 
external field <j>(r), which we will initially assume is nonzero only in a local region, e.g., the potential 
arising when a LJ or hard core solute is fixed at the origin. We describe the system using a grand 
ensemble with fixed chemical potential fi B , which determines p B , the uniform bulk fluid density far 
from the perturbation where 0(r) = 0. The LJ pair potential w(r) = uo{r) +u\(r) is separated into 
rapidly and slowly varying parts associated with the intermolecular forces (67-69) so that all the 
harshly repulsive forces arise from uq and all the attractive forces from u\. Thus uq(t) = w(r) + e 
for r < ro, where r = 2 1 / f V is the distance to the potential minimum where the LJ force changes 
sign, and is zero otherwise, while u\{r) = — e for r < tq and equals w(r) otherwise. Here a and 
e are the usual length and energy parameters in the LJ potential. With this separation u\(r) is 
relatively slowly varying and smooth, with a continuous derivative even at r$. We will make use of 
these features in the theory described below. 

2 Molecular field approximation 

2. 1 Structure of the nonuniform reference fluid 

We now begin our derivation of the VDW interface equation. The fundamental approximation in 
this interpretation of the VDW theory is the introduction of an effective single particle potential 
or "molecular field" that describes the locally averaged effects of the attractive interactions (8) 
in the nonuniform LJ fluid. Since the attractive interactions are relatively slowly varying, such 
an averaged treatment seems physically reasonable. The structure of the resulting "reference" 
or "mimic" system, where the attractive intermolecular interactions have been replaced by the 
effective single particle potential, is supposed to accurately approximate that of the original system. 
Thus the structure of the nonuniform LJ fluid is assumed to be given by that of the simpler 
nonuniform reference fluid (1-3), with only repulsive intermolecular pair interactions uo(r) (equal 
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to the LJ repulsions) and a chemical potential /j,q corresponding to the same bulk density p B but 
in a different renormalized or effective reference field (ERF) 0#(r). Because uq is harshly repulsive, 
many properties of the reference fluid can be accurately approximated by those of a fluid of hard 
spheres with a diameter chosen by the usual "blip function" expansion (10, 66), as described in 
detail in (5, 6). While this approximation is not essential, it is numerically very convenient, and 
for most purposes we will treat the reference system as a hard sphere fluid in the presence of the 
ERF. 

Before we discuss the specific molecular field equation ||] usually used to determine the ERF, 
let us consider some general physical consequences arising from the use of any ERF to describe 
the effects of attractive interactions that will be important in what follows. Since the goal is to 
produce structure in the reference fluid approximating that of the full fluid to the extent possible, 
it is natural to choose 4>R,( r ) m principle so that the local (singlet) densities (70) at every point r 
in the two fluids are equal: 

p o (r;[M/^) = p(r;[0],/i B ). (1) 

Of course this density is not known in advance, so in practice we will make approximate choices 
for (j>R motivated by molecular field ideas. Here the subscript denotes the reference fluid, the 
absence of a subscript the LJ fluid, and the notation [(f)] indicates that the correlation functions 
are functionals of the external field. Unless we want to emphasize this point, we will suppress this 
functional dependence, e.g., writing Equation [l] as po(v) =p(r). 

But the reference fluid can also describe certain features of higher order correlation functions. We 
expect that when 4>r is chosen so that Equation |l| holds, this will produce similar local environments 
for the (identical) repulsive cores in the two fluids, which at high density will mainly determine 
higher order density correlations through excluded volume effects. Thus, when Equation |j] is 
satisfied, it seems plausible that pair correlations are also approximately equal (1, 2), so that 

^(n.^^^Cn.ra). (2) 

In the dense uniform LJ fluid it is well known that correlations are dominated by the repulsive 
forces, and this near equality lies behind the success of perturbation theories of liquids (10, 68, 
69). (However, this structural assumption is rigorously true only in the artificial "Kac limit" where 
the attractive interactions are infinitely weak and long-ranged (64, 71).) The most general use of 
this idea asserts that all structural effects of attractive forces in the nonuniform LJ fluid can be 
approximately described in terms of the structure of the reference fluid in the appropriately chosen 
ERF. 

There are many advantages in using the simpler reference system to define structure. In particu- 
lar, the uniform reference fluid is well defined for all densities from dilute vapor to dense fluid and 
there are no conceptual problems that arise from densities in the two phase region, as would be the 
case in the traditional theories requiring such densities in the original LJ fluid. See, e.g., Section |3| 
below. Moreover, as we will discuss in detail later, there exist simple and accurate theories for the 
reference fluid structure based on linear response theory. 

However there are inherent limitations in such an approximation, most notably in describing 
long wavelength correlations such as those found near the critical point. Some other shortcomings 
of the molecular field approximation when applied to realistic fluid models are easily understood. 
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For example, the effects of long wavelength capillary wave fluctuations (72) of the free liquid-vapor 
interface clearly cannot be described using such reference system correlation functions however 4>r 
is chosen. 

2.2 Simple molecular field equation 

We now turn to the choice of 4>r. In our interpretation, the classical VDW theory uses the sim- 
ple molecular field (MF) approximation for the ERF, given in Equation ||| below. This is just a 
transcription of the usual molecular field equation for the Ising model to a continuum fluid with 
attractive interactions u\{r) and can be arrived at in a number of different ways (8, 44). Particu- 
larly relevant for our purposes here is the derivation discussed in detail in (1-3, 7) that starts from 
the balance of forces in the reference and LJ fluids as described by the exact Yvon-Born- Green 
(YBG) hierarchy (10) and uses Equations |l] and || to arrive at the MF equation by an approximate 
integration. This can be looked on as a modern version of the closely related calculation VDW 
originally carried out (9). The final result is well known: 

<M r i) = <H r i) + J dr 2 po(r 2 ; [4>r], p B ) U!(r 12 ) + 2ap B , (3) 

where 

a = -- J dr 2 u 1 (r 12 ) (4) 

corresponds to the attractive interaction parameter a in the uniform fluid VDW equation, as 
discussed below. The last term in Equation || represents a constant of integration in the derivation 
in (1-3, 7) and is chosen so that vanishes far from a localized perturbation where the density 
becomes equal to p B . 

Because of the integration over the slowly varying attractive component of the intermolecular 
potential ui(ri2), the second term on the right in Equation || is smooth and relatively slowly varying 
even when po itself has discontinuities and oscillations, as could arise from a eft with a hard core. 
Thus the ERF </>ij(ri) is quite generally comprised of the original field (f>(ri) and an additional 
slowly varying term 

4> s (ri) = J dr 2 p (r 2 ) ui(r u ) + 2ap B (5) 

that takes account of spatial variations of the attractive interactions, i.e., the unbalanced attractive 
forces in the nonuniform LJ fluid. 

Some thermodynamic implications of the MF approximation can be seen when we consider the 
MF equation |3| in the limit of a constant field. In the grand ensemble, this is equivalent to a shift 
of the chemical potential, producing a shift in the uniform density, as discussed in more detail 
in the next Section. Thus we can use Equation |3| to relate chemical potentials in the reference 
and LJ systems. Let p(p) and po(p) denote the chemical potential as a function of density in the 
uniform LJ system and reference system respectively. (These also depend on the temperature, but 
we are usually interested in density variations along particular isotherms, so we will not indicate the 
temperature dependence explicitly.) Then (j) p = pP — p(p) is the exact value of that constant field 
such that the density in the LJ fluid changes from p B to p, and <j> p R = p B — po(p) is the analogous 
field in the reference fluid producing the same density change. Using Equation |], the MF equation 
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H gives the MF approximation relating cf> p and (j) R , or equivalently, the MF approximation relating 
p{p) and po(p)- This can be written in the familiar VDW form (64): 

Kp) = W>(p) ~ 2a P- ( 6 ) 

In the limit of a uniform system, Equation || describes all effects of attractive interactions in terms 
of the constant parameter a. Indeed the theory then reduces to the generalized uniform fluid VDW 
theory of Longuet-Higgins and Widom (64, 65), where one combines an accurate description of the 
uniform (hard sphere) reference system with the simple treatment of the attractive interactions in 
terms of a constant background potential a. In the MF approximation p{p) is defined by the right 
side of Equation |6| and has meaning even for values of p in the two phase region. 

In general, to determine 4>R a self-consistent solution of Equation || is required, since <pR appears 
explicitly on the left side and implicitly on the right side through po(r; [(/>r],^q)- In principle, since 
a unique density po(r; [4>r],Pq) is associated with a given external field (pn through the partition 
function, Equation ^ is self-contained and hence self-consistent values for both 4>r and po can be 
found, by iteration, for example. Such solutions were found in (1, 2), using computer simulations 
to accurately determine the associated density po(r; [4>r],Po) for a variety of external fields. 

In practice, one must make additional approximations beyond the molecular field assumption 
to obtain po( r i [<Ar]j Po) m an accurate and computationally practical way. It is here that the 
main limitation of the classical VDW theory arises. In our derivation the classical VDW theory 
results when a second approximation, appropriate for a slowly varying density field (8), is used to 
determine the density po( r ; [4>r]jPo) induced by 4>r. This approximation takes account only of the 
local value of the field through a shifted chemical potential and can be used for a general 4>R( r )- 
When 4>r(t) is very slowly varying it is exact. But in more general cases it gives inaccurate results, 
spoiling the predictions of the VDW theory. To see how this comes about, we first define this local 
hydrostatic density response for a general field, and then show how its use transforms Equation [3| 
into the VDW interface equation as it is usually presented. 

3 Hydrostatic density 

3.1 Local response to general field 

Consider first a given general external field 4>R- Since the chemical potential acts like a constant 
field in the grand partition function, the associated density /?o(r; [4>r],Po) = Po( r ) is a functional 
of 4>R and a function of p^ and depends only on the difference between these quantities. Thus for 
any fixed position ri we can define a shifted chemical potential /Xq 1 = p^ — 4>r(ti) and shifted field 
$>jFj (r) = 4>R(r) — 4>r(ti), whose parametric dependence on ri is denoted by a superscript, and we 
have for all r the exact relation po{r; [(^r],Po) = Po( r 'i [4>r]i Po 1 )- 

By construction the shifted field 0#(r) vanishes at r = ri. If <j>R is very slowly varying, then it 
remains very small for r near rx- In that case, to determine the density at ri we can approximate 
po (r i ; [4> r ^ ] , Pq 1 ) by po (ri ; [0] , //q 1 ) = Pq 1 , the density of the uniform fluid (in zero field) at the shifted 
chemical potential ^Ug 1 . This defines p^ 1 , the hydrostatic density response to the field 4>R- Note that 
Pg 1 depends only on the local value of the field <j)R at ri through its dependence on /ig 1 . When the 
field varies slowly enough the hydrostatic approximation, where |0o( r i) at each ri is replaced by the 
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corresponding uniform fluid density Pq 1 , is very accurate (5). 

One can equivalently define the hydrostatic density /Jq 1 using p,o(p), the chemical potential of the 
uniform reference fluid as a function of density; it is defined by 



This equation plays a central role in all that follows. The hydrostatic density is easy to determine 
for a general 4>R- I n particular if <pR has a hard core with = oo in a certain region of space, then 
the hydrostatic density p^ 1 correctly vanishes in that same region, corresponding to the vanishing 
density of the uniform fluid at the chemical potential p r Q l = — oo. However, because of the strictly 
local response to the field, any nonlocal excluded volume correlations induced by the hard core 
potential outside the hard core region are not properly described by the hydrostatic approximation. 
Still, it remains well defined even in this limit whereas approximations based on using properties of 
a uniform fluid evaluated at the local density (which can easily exceed close packing) break down 
completely (43, 44). 

3.2 Local response to ERF 

To obtain the VDW and MVDW theories we now determine the local hydrostatic response to the 
ERF in Equation ||. This arises after Equation || for <j)R(v\) is substituted into Equation ^: 



As written, Equation ||] just defines the hydrostatic density p^ 1 in terms of the local value of the 
ERF at ri, which itself involves an integral over the full density /9o( r 2; [<Ar]> Pq) at all other points 
T2- To explicitly determine p^ 1 we need to specify Apo(r2) = Po( r 2) — Po 2 , as is clear when Equation 
S is exactly rewritten in the form: 



As we now show, different approximations for Apo( r 2) immediately give the classical VDW and 
the MVDW interface equations. 

4 Hydrostatic approximation and the classical VDW theory 

The classical VDW interface equation arises when one assumes that the full density response to 
the ERF /9o( r 2) for all r-2 is accurately approximated by p^ 2 , the local hydrostatic response. This 
hydrostatic approximation is consistent when the ERF and the induced density are slowly enough 
varying. Thus in this derivation of the classical VDW theory the last term in Equation ^ is ignored 
and Equation || or || is approximated by an integral equation involving only the hydrostatic density: 



Poipl 1 ) = Po ~ <fo(ri). 



(7) 





(9) 




(10) 
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along with the assumption that po(r) = p^ for all r. Here we have also used Equation |6]to replace 
Po(p) by p(p), with the understanding that the latter is really defined by the right side of Equation 
^. Consistent with the assumption of a slowly varying profile and the fact that u\ is reasonably 
short-ranged, Pq 2 in the last term on the right is often expanded to second order in a Taylor series 
about /Jq 1 , yielding in the simple case of a liquid- vapor interface with planar symmetry and (f> = 
a differential equation for the interface profile p§: 

p{p z Q ) = p B + md 2 pl/dz 2 (11) 

where 

m = / dr r 2 u\(r). (12) 
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While there is a strictly local response to the value of the ERF at ri in Equation [ll] yielding 
Pq 1 , the ERF itself at ri depends nonlocally on p r ^ through the integration over the attractive 
interactions. In the classical VDW theory this provides the only source of nonlocality. 



Equations 10 and [II] are completely equivalent to the VDW theory for the density profile as it is 
usually presented (8). In our derivation the theory describes hydrostatic densities in the reference 
system, and p{p) is also defined in terms of reference system quantities given on the right side of 
Equation ||. This provides a simple and consistent interpretation that avoids all problems associated 
with densities in the two phase region of the LJ fluid. 

However, in view of Equation p], one can replace p 1 ^ 1 by p Tl in Equations [h] and 11 and formally 
eliminate all explicit mention of the reference system. This is the way the standard theory is 



usually presented. This formal rewriting of Equations 10 and 11 seems to require only properties 
of the original LJ system and could be useful when applying the theory to other fluids where 
the appropriate reference system is not so apparent. Work for water (4) along these lines will be 
reviewed later. However, this obscures some of the physical underpinnings of the theory in terms of 
the basic MF approximation and it is not clear from the form of these equations how they should 
be corrected in cases where the density profile is more rapidly varying. 

Of course, one might hope that there could exist a more general formulation of the theory where 
one does not need the MF approximation at all. Despite some formal results from density functional 
theory, which we will discuss later, we believe in practice this is not likely to be the case. Long- 
standing conceptual problems (8) arise in interpreting what is meant by p-(p ri ) and p ri itself in 



Equations 1C and 11 for density values in the two phase region of the uniform LJ fluid unless they 
are defined by using a MF approximation either explicitly as in the right hand side of Equation |(| 
or implicitly by assuming some kind of analytic continuation of values from the stable phases. See 
Section || below for further discussion. 

Thus to improve the classical VDW theory, we return to Equations |8| and ||. Here the essential 
physics involving the interplay between the MF approximation for the ERF and the hydrostatic 
density response is clear. If we can understand in detail how to improve the VDW theory for the 
LJ system, we may gain insights that could apply more generally. 



5 Correcting the hydrostatic approximation 
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5.1 Optimized linear response and the HLR equation 



As emphasized in Section 3.2, the VDW theory determines the local density response pjj 1 to the 
ERF 0/j(ri) in Equations || and |9|. Problems arise in the classical VDW theory from the second 
approximation of ignoring the density difference Apo(r 2 ) = Po(r 2 ) — Pq 2 , produced by definition 
from nonzero values of the shifted field 0^ (r) away from r 2 . 

But one can calculate Apo in a very simple way by using linear response theory, which exactly 
relates small changes in the density and field. This approach is clearly correct when Apo and cp^ 
are uniformly small, and in this sense is analogous to a gradient correction to the local density 
approximation in DFT. However, in contrast to the gradient correction, we will see that there are 
good physical reasons to believe that our linear response treatment could remain accurate even for 
large perturbations (3, 5, 6). This will allow us to develop a new and generally accurate theory for 
the nonuniform reference fluid in an arbitrary external field and also will help us correct one of the 
major shortcomings of the classical VDW theory for fluids with attractive interactions. 

We start from the general linear response equation (10) for a nonuniform system in a general field 
4>r, with chemical potential Pq , inverse temperature (3 = (ksT)~ l and associated density po(r): 

- Wfl(ri) = J d^Xo 1 ^!,^, [po])fyo( r 2), ( 13 ) 

which relates small perturbations in the density and potential through the (inverse) linear response 
function 

Xo 1 (ri,r 2 ; [po]) = <5(ri-r 2 )/po(ri)-c (ri,r 2 ; [f>o]). (14) 

Here Co is the direct correlation function of the system with density po(r). In most cases we will 
consider perturbations about a uniform system, soco will take the simple form co(ri 2; p). Since we 
want to focus on effects of the perturbing field, we have used the inverse form of linear response 



theory, where the field appears explicitly only on the left hand side of Equation |13> evaluated at 
ri. 



Could Equation |13| also be used to determine the finite density response Apo to a large field 
perturbation A</>#? Such a linear relation between a (possibly infinite) external field perturbation 
on the left hand side and the finite induced density change on the right must certainly fail for values 



of ri where the field is very large. Conversely, Equation 13 should be most accurate for those values 
of ri where the field is small — in particular where the field vanishes — and then through the 
integration over all r 2 it relates density changes in the region where the field vanishes to density 
changes in the regions where the field is nonzero (3, 5). 

This optimal condition for the validity of linear response theory holds true automatically when 
we use Equation 13 to determine the change at each ri from the uniform fluid with density p^ 1 
induced by the shifted field. Thus we set Xo 1 = Xq 1 ( r i2', Po 1 ) in Equation and take 54>r = </> T ^ 
and <5po(r 2 ) = Po(*"2) — Po 1 - Since <^j(r) is zero a ^ r i by construction, the left side of Equation 13 
vanishes, and we have 

= Jdr 2 XoHr 12 ;pZ)[p (r 2 )-pX], (15) 



which can be rewritten using Equation 14 as 

po(ri) = p'o 1 +P1 1 Jdv 2 c G {r 12 -p^)[p (v 2 )-p^\. (16) 
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This is our final result, which we refer to as the hydrostatic linear response (HLR) equation. A 
more formal derivation involving a coupling parameter integration, and yielding another related 
equation for po( r i) is given in (5). Note that the field appears only implicitly through its local 
effect on pg 1 . The HLR equation is useful by itself as a way to determine the density change 
induced by a known external field <Pr(t). It builds on properties of the uniform reference system, 
requiring in particular co (ri2 ; p) and po (p) . Quantitatively accurate and computationally convenient 
approximations for these functions are known from the GMSA theory of Waisman (73, 74); in 
many cases the simpler Percus-Yevick (PY) approximation (10, 62) discussed below gives sufficient 



accuracy. Given these functions and a known external field <^(r), Equation 16 is a linear integral 
equation that can be solved self-consistently to determine the induced density po( r i) for all ri. We 
will first discuss some general properties of the HLR equation, and examine its accuracy in several 
test cases where the external field ^(ri) is known. Then in Section [O] we will discuss its use in 
correcting the classical VDW theory where 0_r(i"i) is determined self-consistently. 

5.2 Properties of the HLR equation for fixed (fin 



Equation 16 relates the density /fo( r i) at each ri to an integral involving the density po( r 2) at all 
other points and a locally optimal uniform fluid kernel co(ri2; Pq 1 ) that depends implicitly on ri 
through p^ 1 . This ri dependence is the most important new feature of the HLR equation and it 
represents the main reason for improved results over conventional methods, which generally use 
only co(ri2', Pq) along with various nonlinear closures that try to directly relate the field and the 
density in regions where the field is nonzero (10). See, e.g., Equation |17| below. 

5.2.1 Hard sphere solute and the PY equation 

Consider in particular the important test case where (f>R represents the potential of a hard sphere 
solute fixed at the origin. The local hydrostatic response p^ 1 is simply a step function in this case, 
equal to zero for ri inside the core region and to p$ outside. In reality at high density there are 
very large nonlocal oscillatory excluded volume correlations in the true po(ri). As we will see, the 
HLR equation gives results in this special case equivalent to the generally accurate Percus-Yevick 
(PY) approximation. In addition, we can use results from recent computer simulations (75) and 
related work on the Gaussian field model (76) to provide us with further insights into why this 
simple linear response treatment can remain surprisingly accurate even for hard cores. 

We can connect the HLR equation with conventional integral equation theory leading to the PY 
equation by noting that the usual solute-solvent direct correlation function Cs(ri) can be defined 
by 

C s (ri) = J dr 2 xv\rx2;p$)[pv{*2)-p$], (17) 

where exact values for all functions are used on the right side. Thus —Cs(ri)/f3 is that function 
that must replace the perturbing solute-solvent potential to give exact results when the full induced 
density change from p$ is used in the linear response equation |l3[ (Using Equation 14 , we see that 



Equation |17] is equivalent to the standard solute-solvent Ornstein-Zernike equation, which is the 
usual way (10) of defining Cs( r l)-) 
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Similarly, in the derivation of the HLR equation we can replace the zero on the left side of 
Equation 15 by a new function Cs(^i) that exactly satisfies: 

C s (v 1 ) = Jdr 2Xo \r 12 ;p^)[ Po (r 2 )-p^}. (18) 
Note that the local value of the perturbing field <j) r ^ associated with Cs (ri) always vanishes at 



ri, unlike the case for Cs(?i) above. The HLR equation 16 follows from Equation [L8] by setting 
/°o 1( -^s( r i) ec t ua l to zero everywhere, and to the extent that the HLR equation is accurate, we expect 
p^Csivi) to be generally very small. In particular the HLR equation automatically satisfies the 
core condition po = inside the hard core region where both p^ 1 and /5o( r i) vanish. Corrections to 
the HLR equation arise from nonzero values of p^Csi^i) outside the core. 

If we compare Equations [17], with the core condition imposed, and 15 for ri outside the core 
region where pg 1 = p$ , we see they are identical. Thus Cs(ri) outside the core exactly equals 
Cs (ri), the "tail" of the hard sphere solute-solvent direct correlation. This is generally believed 
to be small and in the very accurate GMSA equation it is approximated by a rapidly decaying 
Yukawa function (73, 74). In the HLR approximation this tail is set equal to zero. Thus the HLR 
equation predicts that the hard sphere solute-solvent direct correlation vanishes identically outside 
the core. As is well known, this is equivalent to the PY closure for hard core systems (10). A 
GMSA treatment could be used in Equation ^ to correct the HLR results if more accuracy is 
needed. A self-consistent application of these ideas when the solute is the same size as the solvent, 
requiring that Cs = Co = outside the core, yields the standard PY equation for hard spheres as 
interpreted by Stell (77), with co(r) equal to zero outside the core and po{r) = p$ g®{r) equal to 
zero inside. Similarly, the HLR equation reduces to the PY wall-particle equation when the radius 
of the solute tends to infinity and the PY Co is used. 

5.2.2 Computer simulations and Gaussian fluctuations 

Further insight into the surprising accuracy of linear response theory for uniform hard sphere fluids 
comes from recent computer simulations by Crooks and Chandler (75), following related work 
on water (17). They have shown that even large spontaneous density fluctuations in a uniform 
hard sphere fluid can be accurately described using the same Gaussian probability distribution 
that describes small fluctuations. Small density changes induced by a small perturbing field are 
described by the linear response equation [l3|. By the fluctuation-dissipation theorem, the same 
linear response function Xo~ 1 ( r i2j/ 9 ) controls the spontaneous Gaussian density fluctuations in the 
uniform fluid (10, 76). In particular, they considered the probability distribution for finding N 
hard spheres in a spherical volume of the fluid and found that even cavity or void formation with 
N = was reasonably well described by the Gaussian theory. Since a cavity influences the rest of 
the fluid in exactly the same way as a hard sphere solute of the same size, the density response 
of the hard sphere fluid to such a solute (or imposed region of zero density) should indeed be well 



described using linear response theory with the uniform fluid response function, as in Equation 15 
Thus a key feature exploited in Equation [l6| is that density fluctuations in the reference fluid are 
to a remarkable extent Gaussian. The simulation results pertain to fluctuations in a uniform fluid. 
Similarly, to describe effects of a general external field, the HLR equation considers perturbations at 
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each ri to the uniform hydrostatic fluid density Pq 1 induced by the shifted field (fy , which vanishes 
at least locally at r\. 



5.2.3 Other rapidly varying model potentials 

Thus we see that the HLR equation ^ satisfies the following limits, i) It is exact when 4>n(r) 
is very slowly varying and exactly describes to linear order small corrections to the hydrostatic 
approximation ii) It is exact for any 0(r) at low enough density, where there is a local relation 
between the potential and induced density, iii) For a field c/>(r) from a general hard core solute, 
Equation |0] reduces to the PY approximation, as discussed above. To give more indications of the 
accuracy of Equation 1€ , we review solutions (5) for some model potentials designed to show both 
the strengths and the weaknesses of the present methods and compare with computer simulations. 

Figure 1 shows the correlation functions go( r ) = Po{ r )/Po f° r a hard sphere system at a mod- 
erately high bulk density p$ = 0.49 induced by the deep attractive spherical parabolic model 
potentials shown in the inset. (Reduced units, with distances measured in units of the hard sphere 
diameter are used.) The HLR equation reproduces the increased density inside the well, and the 
nonlocal oscillatory excluded volume correlations, which show a local density minimum at the 
center of the well due to packing effects. 

Since the external field enters Equation [T6| only locally through its effect on p 1 ^ 1 , it is also easy 
to use it for the inverse problem of determining the field associated with a given density profile. 
As an example, the crosses in the inset gives the potentials predicted by Equation 16 given the 
simulation data for g(r). 

While these results are qualitatively very satisfactory for the most part, there are some problems. 
Note in particular the slightly negative density at the center of the well; a positive density is not 
guaranteed in this linear theory. (The related HM equation derived in (5) always yields a non 
negative density and does slightly better here but it performs less well in the hard wall limit.) One 
would expect methods based on expanding about the uniform hydrostatic fluid to be least accurate 
for potentials with very steep gradients, and poor results at high density were found for a repulsive 
planar triangular barrier potential that rose from to lOfc^T over a distance of one hard sphere 
diameter. Results actually improve as the barrier height increases and the potential approaches a 
hard wall potential: the theory does better for hard cores because there is no region of space where 
there is a large gradient in the external field while at the same time the local density is nonzero. 
A GMSA like treatment based on Equation 18 may improve matters here. But large repulsive 
potentials with very steep gradients are better treated by "blip function" methods (66) or other 
expansions about a hard core system. 



6 Correcting the molecular field approximation 

We now return to fluids with attractive interactions as described using a self-consistently determined 
ERF. Thinking about the hydrostatic limit also can help us correct some quantitatively inadequate 
features of the simple MF approximation ^ for the ERF (7). In the limit of a uniform system this 
describes all effects of attractive interactions in terms of the constant VDW parameter a. While this 
very simple approximation captures much essential physics and gives a qualitative description of the 
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uniform fluid thermodynamic properties it certainly is not quantitatively accurate. For example, 
when used to describe a slowly varying liquid-vapor interface, it will predict shifted (molecular 
field) values for the densities of the coexisting bulk liquid and vapor phases. The main problem 
with the theory is not so much its description of the local density gradients, but its predictions for 
the thermodynamic properties of the coexisting bulk phases (60). 

To achieve quantitative agreement with bulk thermodynamic properties one can replace the 
constant a by a function a that depends (hopefully weakly, to the extent the classical VDW 
theory is reasonably accurate) on temperature and density. This has been suggested many times 
in the literature, usually in the context of perturbation theories for uniform fluids (10), where 
approximate expressions giving such functions have been derived (78). We review here a simpler 
and more empirical approach (7) that can be used to determine a if an accurate analytic equation of 
state for the bulk fluid is known, and we will use this in a simple modification of the MF expression 
for the ERF to insure that exact thermodynamic results (consistent with the given uniform fluid 
equation of state) are found in the hydrostatic limit of a very slowly varying ERF. 

To that end, let us consider again the MF equation in the limit of a constant field. Instead of 
using the MF approximation for /i(p) as in Equation || we assume that p(p) is known from an 
accurate bulk equation of state. In particular, we determine p{p) from the 33-parameter equation 
of state for the LJ fluid (79) given by Johnson, et al. This provides a very good global description 
of the stable liquid and vapor phases in the LJ fluid and provides a smooth interpolation in between 
by using analytic fitting functions. Thus it naturally produces a modified "van der Waals loop" in 
the two phase region and seems quite appropriate for our use here in improving the simplest MF 
description of the uniform fluid. 

Using known properties of the hard sphere fluid, we also have an essentially exact expression for 
Po(p)- Using this we can relate p(p) and po(p) in analogy to Equation [6|: 

p{p) = p (p) - 2a(p)p (19) 

but with the constant a replaced by a function a(p) of density and temperature chosen so that 



Equation |19| holds. Since even the simplest molecular field approximation is qualitatively accurate 
we expect that the ratio a(p)/a will be of order unity and rather weakly dependent on density and 
temperature. 

We indeed found that the constant a was a good overall compromise (7). However the true a(p) 
fit to the equation of state exhibits variations of up to about fifteen per cent as a function of density 
and temperature, illustrating the need for an accurate equation of state for quantitative accuracy. 

Because of the strictly local response, these results for a constant field can be used to determine 
exact results in the hydrostatic limit of very slowly varying fields. We want to modify Equation 
^ so that in the hydrostatic limit it will reproduce these exact values, while still giving reasonable 
results for more rapidly varying fields. 

There is no unique way to do this, but the following simple prescription seems very natural, and 
gives our final result, which we will call the modified molecular field (MMF) approximation for the 
ERF: 

fl ( ri )-0(n) = ^^ / dr 2 p (r 2 ;[ ( PR}^o)Mn2)+2a(p B )p B . (20) 
Thus the molecular field integral in Equation || is multiplied by a factor a(p 1 [ j 1 )/a of order unity 
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that depends on ri through the dependence of the hydrostatic density pg 1 on the local value of the 
field (j)ji(ri), and the constant of integration 2ap B is replaced by the appropriate limiting value of 
the modified integral. Note that the hydrostatic density Pq 1 remains smooth and relatively slowly 
varying outside the core even when 0#(ri) contains a hard core. The nonlocal oscillatory excluded 
volume correlations that can exist in the full density /Oo( r i) do not appear in p^ 1 because of the 
strictly local response to 4>r. Results using Equation ^0] will be discussed below. 



7 MVDW theory of nonuniform fluids with attractive interactions 

7.1 Two step method 

The ERF must be calculated self-consistently for fluids with attractive interactions. Our new 
MVDW theory corrects the classical VDW theory by using Equation ^ to accurately determine 
po( r 2) in Equation ||. Thus the MVDW theory requires the simultaneous solution of two equations: 
Equation || (or the closely related equation that arises if the Equation 2C is used for the ERF) , and 
Equation [u]. Equation || determines p 1 ^ 1 , the local hydrostatic response to the ERF as in the VDW 
theory described above, and Equation determines the full nonlocal response po( r i)- (Similarly, 
the VDW theory can be viewed as replacing Equation 16 by the hydrostatic approximation po( r i) = 
Pi 1 -) 

One can think about solving these equations by a two step iterative method (3). For a given 
approximation to the ERF, determine in a first step the associated smooth hydrostatic density p 1 ^ 1 
from Equation ||. Then in a second step take account of nonlocal and usually oscillatory corrections 
to this profile, generally induced by rapidly varying features in the external field (j), using a locally 



optimal application of linear response theory, Equation 16. This new density po( r i) is then used 
to compute a new approximation for the ERF, and the two steps are iterated to self consistency. 
This process is easy to implement numerically and rapid convergence has been found in most cases 
we have examined (3, 6, 7). 

We noted in Section that the simple MF ERF can quite generally be written as 4>r = (j) + (f> s , 
where (f) s in Equation || takes account of the unbalanced attractive interactions and is smooth and 
relatively slowly varying. We are often interested in cases where <j) is a hard core potential. (If <j) 
also has a slowly varying part, say from weak attractive interactions between a solute and solvent, 
the latter should be added (6) to <p s in the discussion that follows.) In this case one can implement 
the two step process in a slightly different way, which is physically suggestive and was in fact how 
we first thought about the problem (3, 4, 6). 

Let us consider in the first step the hydrostatic response p r / to the slowly varying part (ft s of the 
ERF alone. Since there is a strictly local response to the ERF, p^ 1 differs from p^ 1 for the full ERF 
only inside the core region, where p 1 ^ 1 vanishes while p* 1 remains continuous and smooth. This 
smoothness allows us to use the gradient approximation in the next to last term of Equation ^ in 
determining p* 1 if we wish. Since <p s describes the unbalanced attractive interactions associated 
with interface formation, we can interpret the smooth p* 1 as an interfacial component of the full 
density response po(r\). 

Then in the second step we take into account the response to the remaining hard core part of 4>r 
and all nonlocal effects. This will cause po(?i) to vanish inside the core and, depending on the extent 
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to which the density p^ 1 near the core has been reduced, can induce nonlocal oscillatory excluded 
volume corrections to p* 1 outside the core, which we again calculate using linear response theory. 
Thus speaking pictorially, the second step takes into account the nonlocal Gaussian fluctuations 
induced by the hard core potential about the smoothly varying interfacial component of the density 
profile. The strongly non-Gaussian component associated with interface formation and arising from 
the unbalanced attractive interactions is taken into account in the first step. The final profile results 
from the self-consistent interplay between the components described in each step. 

Since we are using linear response theory in all cases to correct a hydrostatic response, these 
alternate ways of implementing the two step procedure are completely equivalent. In particular in 
the first interpretation described above we calculate the hydrostatic response Pq 1 to the full ERF 
and correct it for all ri using Equation [l^. The hard core condition po(r) = for all r inside 
the core region follows automatically from perturbing about the local hydrostatic density. In the 
second interpretation we calculate in the first step the hydrostatic response pi 1 to the interfacial 
component 4> s alone and correct it by using linear response only for ri outside the core, while 
imposing /9o( r 2) = in the core region: 







/draXoVia^M^-pI 1 ]- (21) 



Since p r 8 x equals Pq 1 outside the core, the solution po( r ) to Equation 21 with the core condition 



imposed is identical to that given by Equation 15. 



7.2 Hard sphere solute in L J fluid near liquid-vapor coexistence 

These ideas have been applied to the nonuniform LJ fluid in a variety of different situations, 
including fluids near a hard wall, fluids confined in slits and tubes, and the liquid-vapor interface, 
with generally very good results (1-3, 5-7, 80-82). We will review here results (7) for a system 
studied by Huang and Chandler (HC) that combines many of these limits (83), and is physically 
very relevant for our later discussion about hydrophobic interactions. 

HC carried out extensive computer simulations to determine properties of the LJ liquid at a 
state very near the triple point with p B = 0.70 and T = 0.85 as the radius of a hard sphere solute 
is varied, and compared the results to the LCW theory reviewed below. (We use the standard 
LJ reduced units.) By definition the solute centered at the origin interacts with the LJ particles 
through the hard core potential: 

«r,S)-(~- (22) 
[0, r > b. 

The MMF theory discussed above allows us to reduce this problem to that of the reference fluid 
in the presence of the effective field (j)R(r;S) satisfying Equation We have calculated self- 
consistently the ERF (j>n(r; S) and the associated density response po(r; S) of the reference fluid, 
solving Equations ^0] and |0] by iteration. In Figure 2 we compare these results for the density 
profiles in the presence of the hard sphere solutes with S equal to 1.0, 2.0, 3.0 and 4.0 in reduced 
units with the simulation results of the same LJ system by HC. There is very good agreement 
between theory and simulation. 
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Figure 3 shows the corresponding ERF's obtained in these calculations. For small solutes with 
S less than about 0.7, attractive interactions do not give rise to any substantial modification of 
the bare external field, as can be seen from the plot of 4>R,{r; S) for S = 0.5. (Clearly, for S = 
there are no solute induced interactions of any kind and the profile reduces to the constant p B ). 
However, the effects due to unbalanced attractions become important even for S = 1.0, which is 
about the same size as the LJ core, and all larger sizes give rise to a very strong and relatively 
soft repulsion in 4>R(r; S). The corresponding density profiles show pronounced depletion near the 
surface of the solute, characteristic of surface induced drying. 

7.3 Solvation free energy 

Another quantity of great interest is the free energy of the nonuniform system. This is the main 
focus of attention in density functional theory, briefly discussed in Section || below. In contrast, 
the MVDW approach focuses first on the liquid structure. We believe this permits physical insight 
to play a more direct role. However, since we can determine the density response to an arbitrary 
external field, the free energy can be easily calculated from a coupling parameter type integration 
that connects some initial state (e.g., the bulk fluid) whose free energy is known to the final state 
as the field is varied. In the present case there is a very simple route to the free energy of the 
nonuniform LJ system that uses structural features that we know from simulations are accurately 
determined. 

In particular let us consider the change in free energy of the LJ fluid as the range of the external 
field representing the hard core solute is varied from zero to its full extent S. This construction is 
the basis of scaled particle theory (84), and it is well known that the free energy change takes a 
particularly simple form: 

l 

(3Att s = 4vr5 3 J dX X 2 Px (\S + ), (23) 
o 

which requires only the contact value p\(\S + ) of the density profile. This is very accurately given 
by the theory described above. To use this "virial route", we can replace the A-integration by a 
sum and calculate the density profile for several values of A at the fixed bulk chemical potential 

Using Equation ^ we obtain the dependence of solvation free energy on the size of the hard 
sphere solute. The free energy per unit surface area of the solute AQs/^ttS 2 we obtain is shown 
in Figure 4. For small solutes unbalanced attractive interactions do not play an important role 
and the solvation free energy agrees well with a pure hard sphere model, which completely neglects 
attractions by using the bare hard core solute potential, as shown by the dotted line. At the solute 
size of about 0.7 the behavior changes drastically and the reduced free energy rapidly crosses over 
to the practically constant plateau in agreement with the simulation results. The small slope of 
the curves in Figure 4 for large S can be understood by separating the free energy into volume 
(Vs = 47r5 3 /3) and surface (As = 47r£ 2 ) contributions as discussed by HC (83): 

AQ S ~ V s p B + A slS - (24) 
The first term in this expression corresponds to the work required to remove liquid particles from 



Theory of nonuniform liquids 19 

the volume occupied the solute, where p B is the bulk liquid pressure, and is very small for the values 
of S considered here. The second term determines the cost of forming the liquid-solute interface 
and is proportional to the interface tension 75, which is essentially independent of the solute size 
for large solutes. 

8 Hydrophobic interactions in water 

We now discuss an important extension of these ideas to hydrophobic interactions in water (4), as 
described by Lum, Chandler and Weeks (LCW). We first consider water at ambient conditions as 
the radius of a hard sphere solute at the origin is varied. Because this state is very close to the 
liquid-vapor phase boundary, phenomena involving interfaces can be very important. This system 
serves as a simple model of a hydrophobic object in water — the solute does not participate in 
hydrogen bonding and it creates an excluded volume region where the density of water molecules 
vanishes. Weak van der Waals attractions between the solute and solvent do not change the 
qualitative nature of the phenomena we will discuss and can easily be taken into account (6, 85). 

But how can one sensibly apply the theory to water? The local structure of water is certainly 
very different from that of the LJ fluid and anything relying on the detailed properties of a hard 
sphere reference system cannot be trusted. However the two systems do have certain essential 
features in common that can be exploited in a properly generalized MVDW theory. 

First, small fluctuations in liquid water are Gaussian, and computer simulations had earlier 
shown that even relatively large fluctuations leading to the formation of molecular sized cavities 
can be well described using the same Gaussian distribution (17). Indeed Pratt and Chandler (13) 
developed a quantitative theory for the solvation of small apolar molecules using the experimental 
linear response function for water that takes into account the structural and free energy changes 
induced by the excluded volume of the solute. These ideas have been significantly extended and 
placed on a firmer conceptual basis in recent work by Pratt, Hummer and coworkers (17, 19, 20, 
22, 30, 86). See also (76). This suggests that the Gaussian/linear response ideas used in the second 
step of the two step method in the MVDW theory, using response functions appropriate for water, 
could be modified to apply to water. 

Second, just as for the LJ fluid above, the main non-Gaussian feature to be expected in this 
application is associated with interface formation. Long ago, Stillinger (11) gave a qualitative 
and physically very suggestive description of what would be expected as the radius of the hard 
sphere solute is increased. While small solutes should not significantly disturb the hydrogen bond 
network, which can simply go around the solute, in the vicinity of a sufficiently large solute or wall 
the network must be completely disrupted. In the latter case, Stillinger argued, the arrangement 
of molecules near the solute and the interface free energy should resemble that of the liquid-vapor 
interface, which optimally solves the similar problem of going from a complete hydrogen-bond 
network in the liquid to no hydrogen bonds in the dilute vapor. 

Thus Stillinger envisioned a drying transition very like that studied in the last Section for the LJ 
fluid as the solute size is increased. This phenomena is very general. The differences in local struc- 
ture should have little effect on the generic and qualitative physics leading to interface formation. 
Thus it seems plausible that the interfacial component determined in the first step of the MVDW 
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theory could be described at least qualitatively by a MF treatment similar to that used for the LJ 
fluid provided appropriate thermodynamic parameters for water are used. 

We see that key features of both steps of the MVDW theory have some analogues for water. 
The hard sphere reference fluid picture for the LJ fluid was incisive in developing the general ideas 
leading to the MVDW theory. Given that understanding, we may be able to develop an analogous 
approach for water and other fluids that does not rely on the details of the reference fluid, or indeed 
explicitly introduce a reference fluid at all. 

Thus, following LCW, let us examine the essential features of the MVDW theory and see how 
they can be modified to apply to water. As in the last part of Section we will try to describe 
everything formally in terms of the properties of water itself, and not a reference system, though 
molecular field ideas will be introduced to define what is meant in two-phase regions, etc. A crucial 
part of the physics of interface formation in the MVDW theory is the description of the unbalanced 
attractive interactions in Equation |5[ In the LJ fluid this can be rewritten and reinterpreted in terms 
of an averaged or coarse grained density p(ri) with a normalized "weighting function" proportional 
to Uf. 

- 2ap(ri) = J dr 2 p(r 2 )u 1 (r 12 ). (25) 

LCW argue that the unbalanced attractive interactions in water can be described by a similar coarse 
graining of the water density, with the coarse graining carried out over the appropriate range A of 
the attractive interactions in water. 



Now consider the first step of the two step method, as described in the last part of Section 7.1 



where the smooth interfacial component p^ 1 associated with the field <p s in Equation |5| is determined 



from Equation |9|. Using the notation of Equation |25| and expanding the next to last term, Equation 
H becomes: 

MPs 1 ) = M B + mV 2 p r / + 2a[p( ri ) - p*/). (26) 

With appropriate change of notation this is exactly Equation 5 of LCW. However LCW did not con- 
sistently interpret p* 1 (= n s (r\) in LCW) as the hydrostatic density and some later simplifications 
that arise from this were not exploited. Since the slowly varying and generic interfacial component 
p* 1 in Equation ^ should be essentially independent of local structure, LCW used a simple VDW 
form for p(p), which automatically interpolates in the two phase region, but with VDW parameters 
a and b chosen to reproduce the density and compressibility of liquid water at phase coexistence 
and T = 298K. The parameter m was fit to the surface tension of water and the VDW relationship 
between a and m determined the coarse-graining scale A for /5(r), which LCW carried out using a 
simple Gaussian weight. While the details of this fitting procedure are rather arbitrary, and the 
VDW equation is probably inadequate to describe the quantitative relation between the energy 
density and surface tension in water (63), LCW showed that small variations in these parameter 
values did not change the qualitative picture that emerged. 



Now turn to the second linear response step, described for the LJ fluid by Equation 21. What 
is the analogue of the LJ reference Xo 1 ( r i2i p) f° r water? In view of Equation ^, one can formally 
consider fluctuations in the full fluid instead. For p in the stable liquid phase, small fluctuations are 
Gaussian and the uniform fluid x (?"i2i p) can be used directly. However when interfaces form, we 
need an approximation for \~ l in our MF treatment that remains well defined for all values of p in 
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the two-phase region as well as in the vapor phase, as is the case for the reference Xo 1 f° r the LJ 
fluid. In effect LCW devised an interpolation scheme for a x _1 for water that has these properties. 



LCW essentially considered an alternate but equivalent version of Equation 21: 

Pin)-?. 1 = jdr 2X {r 12 -p^)Cs{r2), (27) 

which involves the standard linear response function x( r i2! p) = pd( r i — ? 2) + p 2 h(r 12; p), the inverse 
to x 0*121 p)- Here h+1 is the radial distribution function for water and Cs is a generalized solute- 
solvent direct correlation function, similar to that in Equation |l7]. This is nonzero only inside the 
core (PY closure), and is completely determined by the imposed requirement that p(r) equals zero 
for all r in the core region. They replaced h(r '12; p) by h(r\2\ p B ) in the definition of x> an d used 
experimental values for the latter. This approximation is exact in the bulk liquid phase, reduces 
to the correct ideal gas value at very low density, and smoothly interpolates in between. Equation 
6 in LCW reduces to our Equation ^7] when it is realized that p^ 1 is the hydrostatic density. Since 
p^ 1 is slowly varying, no essential differences should result from using either equation. 

Thus LCW implement in an approximate, but plausible way both steps of the MVDW theory, 
making use of experimental thermodynamic and structural data. They obtained results for the 
solute-solvent density distribution function qualitatively very similar to that of Figure 2 for the 
analogous LJ-hard sphere system. As the radius of the solute increased they found a crossover 
on the biologically relevant length scale of nanometers from "wetting" with peak densities greater 
than the bulk to "drying" with peak densities less than the bulk. (The unbalanced attractive 
interactions cause significant density perturbations in both regimes, so this terminology is somewhat 
misleading.) Weak attractive van der Waals attractions between the solute and water can shift 
inward the position of the interface when partial drying occurs and suppress complete drying, but 
should not have an important effect on the basic interface structure and free energy changes or the 
length scale for the crossover (85). 

LCW also calculated solvation free energies, using a Gaussian approximation that is quantita- 
tively somewhat less accurate than the coupling parameter method discussed above for the LJ 
system, but quite sufficient for qualitative purposes. Again they found behavior very similar to 
Figure 4, with free energies scaling with surface area only on large length scales of order nanome- 
ters. LCW also studied assemblies of extended idealized hydrophobic objects (plates and rods) and 
found that the drying can lead to strong attractions between sufficiently close pairs of such sur- 
faces as the intervening water is expelled. Thus there is a length scale dependence of hydrophobic 
interactions. LCW suggest that such phase transitions could play an important role in aspects of 
protein folding where extended mostly hydrophobic regions approach one another. 

It is beyond the scope of this article (and the expertise of this author!) to assess the validity 
of that last conjecture, given the many complications occurring in nature. Clearly much more 
theoretical and experimental effort is called for on all aspects of the theory. Some interesting recent 
work along these lines can be found in (21-29, 87-93). Our purpose is merely to argue that interface 
formation is a fundamental piece of physics that almost certainly occurs for the idealized models 
discussed herein and that LCW have developed a qualitatively reasonable MF treatment of that 
process based on sound statistical mechanical principles. From that perspective, the scaling of 
hydrophobic solvation energies with exposed surface area with a value close to the surface tension 
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of water can be justified only on large length scales (4, 93). 
9 Density functional theory 

As discussed in the Introduction, the most commonly used theory for the structure and thermo- 
dynamics of nonuniform fluids is density functional theory (DFT). Space limitations permit only a 
few general remarks here focusing on possible relations to the MVDW theory and the advantages 
and disadvantages of each approach. Two standard reviews of DFT are found in (43, 44). 

A main focus of DFT is the intrinsic free energy density functional F([p],[w]), which arises 
from the usual grand canonical free energy Q by a Legendre transform (10), where the functional 
dependence on the chemical potential and external field 0(r) is replaced by a functional dependence 
on the (uniquely associated) singlet density p(r). Our notation emphasizes that F, just like O, 
remains a functional of the intermolecular pair (and any higher order) potentials w. In principle 
this formalism is exact and F contains all the information in f2. An exact hierarchy of (direct) 
correlation functions can be derived by successive functional derivatives of F with respect to p. 

While it must be as hard to determine F exactly as it is Vt, the hope is that practically useful 
approximations may suggest themselves more naturally in this representation, and the sound the- 
oretical basis of DFT may provide a means for systematic corrections. By starting from the free 
energy, certain exact sum rules relating integrals of those correlation functions to the thermody- 
namic properties are automatically and consistently satisfied in any DFT (45, 59, 60); this is not the 
case for structurally based methods such as the MVDW theory. Particularly suggestive is the fact 
that the classical VDW theory can be viewed as a DFT for F where a MF approximation is made 
for the functional dependence on w and a local density approximation is made for the functional 
dependence on p. While this certainly seems promising, it turns out to be rather difficult to do 
significantly better from this starting point. 

Thus, can the DFT formalism help us improve on the MF approximation? Consider the first 
functional derivative of F, which from basic properties of the Legendre transform is easily seen to 
satisfy 

SF([p],[w])/5p(r) = /x-0(r). (28) 

In principle if we knew the "exact" F we could solve this equation, determining the density p(r) 
associated with a given potential 0(r). While most discussions have focused on the density de- 
pendence, the main problem in determining an accurate F from a fundamental point of view is 
its dependence on w. When there are attractive interactions, the exact F must describe critical 
phenomena, capillary waves, and a host of other properties for which we have essentially no idea 
what the true functional dependence on w should be. A treatment incorporating even the simplest 
capillary wave correlations for the liquid-vapor interface produced a density functional very differ- 
ent from conventional ideas (44, 72, 94), and this does not begin to address the range of phenomena 
arising from the general functional dependence on w. So as a practical matter we are essentially 
forced to accept the MF treatment of the attractive interactions, perhaps modified slightly as in 
Section ^ to give a better description of the uniform fluid (60). 

The true F for a system with repulsive forces only must be significantly simpler. Even here the 
dependence on a general repulsive wq can cause problems; any density functional must also explicitly 
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or implicitly approximate the dependence on wo- It is not obvious from a fundamental point of 
view how to modify a functional that gives a good description for say hard spheres, so that it can 
describe much softer repulsions. See (53, 54) for interesting work on this question. Unlike most 
quantum-mechanical applications of DFT, where there is a single coulomb interaction potential, 
we must deal with a range of interactions and hence potentially many different functionals. 

Most workers have quite reasonably focused on the basic hard sphere system, and here some 
significant progress has been made in determining accurate approximations for Fo, particularly in 
the development of fundamental measure theory (50-52, 58). However, the latter is very complicated 
and one might hope for a simpler and perhaps more physically motivated approach. 

Unfortunately, the DFT formalism itself gives few indications of how to proceed in a practical 
manner. The basic problem can be seen when we consider Equation |2^ for a system of hard spheres 
and a general <j)R, i.e., 

SF ([p ], [u ))/dpo(r) =p$- Mr)- (29) 

This equation can be viewed as a generalization of our hydrostatic equation [?], which in principle 
could be used to determine the exact density if a proper Fo could be found. Indeed Equation [7| 
follows from Equation ^ for a very slowly varying profile when a local density approximation is 
made. But how should this be corrected for rapidly varying densities? To go beyond the local 
density approximation, various weighted density approximations have been proposed. However, 
from a fundamental point of view except in special low dimensional cases (58) it is not even clear 
that Fq can be written as a simple functional of weighted densities; certainly the original f?o is 
not a functional of a weighted 0ij(r)! This contrasts with the development of the MVDW theory 
from the VDW theory. There the simplest correction (linear response) is clearly correct for small 
perturbing potentials, and it remains reasonably accurate even in the hard core limit, reducing 
naturally to the PY approximation. Usually in DFT this limit is imposed by hand. 

Since many treatments of DFT have discussed the formal advantages of the method, we have 
mostly focused here on some of the difficulties as we see them. Despite these (perhaps pedantic!) 
objections, there have been impressive successes arising from DFT, and there clearly are close 
connections between some of the ideas of DFT and the MVDW theory. In particular, while we do 
think about molecular fields, the HLR equation 16 explicitly involves only densities. Indeed the 
reason for its success is the removal of any explicit dependence on the field through the expansion 
about the hydrostatic density. However, a deeper connection to DFT has so far escaped us; this is 
an active topic of our current research. 



10 Final Remarks 

We conclude with a few general remarks. The MVDW theory combines in a self-consistent way 
two standard and widely used ingredients: molecular (or mean) field theory and linear response (or 
Gaussian field) theory. The molecular field equation for the ERF takes account of the unbalanced 
attractive forces discussed in the original VDW theory (9); all that is required is to determine 
the molecular field and the induced structure accurately. While the hydrostatic approximation 
used in our interpretation of the VDW theory is not generally accurate, it serves as an optimal 
starting point for corrections based on linear response theory. The resulting theory can handle 
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problems involving fluctuations on a variety of length scales where each ingredient alone would fail. 
Because of its sound physical basis, we are hopeful that the MVDW theory will prove useful in 
many different applications. In particular, assessing and further developing both the physical and 
biophysical implications of the LCW theory in realistic environments seems an important topic 
for future research. One new direction we are thinking about involves fluids with long-ranged 
coulomb forces. Here there are a wide range of phenomena such as charge ordering and double 
layer formation arising from very strong and competing interactions on many length scales. This 
will certainly put our current ideas to a most severe test. 

There are also many basic theoretical questions left open. Can the MVDW theory be understood 
as some kind of weighted DFT? Is it possible to extend these ideas to solid-fluid interfaces? Can one 
go significantly beyond the molecular field picture in describing the effects of attractive interactions 
while still maintaining a tractable theory? The LCW theory, despite its use of properties only of 
water, is definitely a molecular field theory. But the idea of short wavelength Gaussian fluctuations 
related to local structure and long wavelength slowly varying fluctuations related to interface for- 
mation seems more general. A reinterpretation of LCW theory from this perspective is found in 
(63). We close with what is probably the most basic question: Is it is generally correct to imagine 
that fluctuations are essentially Gaussian in most non-critical liquids except when interfaces form? 
A deeper understanding of whether and why this is true is called for. If this physically suggestive 
picture remains valid, then the MVDW theory has captured in a surprisingly simple way much of 
the non-critical physics of the liquid state. 
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Figure 1: Correlation functions for hard spheres in the presence of spherical parabolic potentials 
shown in the inset (solid lines) as given by theory and simulation. The upper curve corresponds to 
the smaller potential and has been displaced upward by one unit. Also shown in the inset (crosses) 
are the potentials predicted by Equation |l^ given the simulation data. 
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Figure 2: Density profiles of the LJ fluid (T = 0.85, p B = 0.70) in the presence of the hard sphere 
solute with 5 = 1.0, 2.0, 3.0 and 4.0. Circles denote simulation results (83). Lines are results of 
the self-consistent approach based on the modified molecular field determined from Equation p0| . 
For ease of viewing, the density profiles for S = 1.0, 2.0 and 3.0 have been shifted vertically by 0.6, 
0.4 and 0.2 units respectively. 
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Figure 4: Dependence of the solvation free energy on the cavity size S. Circles denote results of 
simulations (83) . Lines are obtained from Equation [2^ by using the results of the molecular field 
equation [2(] (solid) , and by neglecting the molecular field (dotted) . 



